---
title: "Main Analysis for Assessing Public Value Failure in Government Adoption of AI"
author: "Daniel Schiff, Kaylyn Jackson Schiff, and Patrick Pierson"
date: "2021"
output: pdf_document
editor_options: 
  chunk_output_type: console
---

#####Setup Chunk#####
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(dplyr)
library(ggpubr)
library(tidyr)
library(forcats)
library(stargazer)
library(sandwich)
library(broom)
library(dotwhisker)
library(rmarkdown)
```

#####Read in the Clean Data####
```{r read in data}
load("Data/ads_clean.RData")
```

####Reproduce Information for Table 1####
```{r table 1}
#Create table
table_1 <- ads %>% group_by(Vignette) %>% summarize(n = length(Vignette))
table_1
```

####Reproduce Information for Table 2####
```{r table 2}
#Perform regressions
out1_adjust <- lm(Feeling ~ Mechanism + Pol + Knowledge + Gender + Age + Race + Educ + Marital_Status + Income + Employment, data=ads)
out2_adjust <- lm(Trust ~ Mechanism + Pol + Knowledge + Gender + Age + Race + Educ + Marital_Status + Income + Employment, data=ads)
out3_adjust <- lm(Quality ~ Mechanism + Pol + Knowledge + Gender + Age + Race + Educ + Marital_Status + Income + Employment, data=ads)
out4_adjust <- lm(Impact ~ Mechanism + Pol + Knowledge + Gender + Age + Race + Educ + Marital_Status + Income + Employment, data=ads)

#Create table
stargazer(out1_adjust, out2_adjust, out3_adjust, out4_adjust,
          se=list(sqrt(diag(vcovHC(out1_adjust,type="HC1"))), 
                  sqrt(diag(vcovHC(out2_adjust,type="HC1"))),
                  sqrt(diag(vcovHC(out3_adjust,type="HC1"))),
                  sqrt(diag(vcovHC(out4_adjust,type="HC1")))),
          omit=c("GenderOther", "Marital", "Employment", "Industry"),
          covariate.labels=c("Bias", "Lack of Transparency", "Lack of Responsiveness",
                             "Liberal", "Conservative", "Low AI Knowledge", "Some AI Knowledge", "High AI Knowledge",
                             "Female", "Millenial", "Gen X", "Boomer", "Silent Gen", "Black", "Other Race", "Some College",
                             "Bachelor's or Higher", "Middle Income", "High Income"),
          title="Public Value Failure Regression Results (Covariate-Adjusted)",
          label="mechanism_table_adjusted",
          notes="With robust SEs, including covariates, and without block coefficients and some covariates printed",
          style="APSR",
          header=F,
          type="latex",
          out="Tables/table_2.tex")
```

####Reproduce Figure 5####
```{r figure 5}
#Create figure
pol_outcomes <- ads %>% group_by(Pol, Vignette) %>% filter(Pol != "Independent") %>% summarize(
  mean_feeling = mean(Feeling), sd_feeling = sd(Feeling),
  mean_trust = mean(Trust), sd_trust = sd(Trust),
  mean_quality = mean(Quality), sd_quality = sd(Quality),
  mean_impact = mean(Impact), sd_impact = sd(Impact))

pol_outcomes <- pol_outcomes %>% pivot_longer(cols = 3:10) %>% 
  mutate(outcome = rep(rep(c("Feeling", "Trust", "Quality", "Impact"), each = 2),8)) %>% 
  ungroup %>% 
  mutate(name = rep(c("mean","sd"),64)) %>% 
  pivot_wider(id_cols = c(1,2,5), names_from = name, values_from = value)
pol_outcomes$outcome <- factor(pol_outcomes$outcome, levels = c("Feeling", "Trust", "Quality", "Impact"))

feeling_pol <- ggplot(pol_outcomes %>% filter(outcome == "Feeling"), 
                      aes(x=fct_rev(Vignette), y=mean, color = Pol, shape = Pol)) + 
  geom_point(size = 1.5) +
  geom_hline(yintercept=0, colour = "grey60", lty=2) + 
  scale_color_manual(values = c("#0015BC","#e9141d")) +
  scale_shape_manual(values = c(17,16) ) +
  coord_flip() +
  xlab("") + 
  ylab("\n") +
  labs(title="Feeling") +
  ylim(-5,5) + 
  theme_bw() + 
  theme(plot.title = element_text(face="bold", hjust = 0.5),
        legend.justification=c(.6, 0),
        legend.background = element_rect(colour="grey80"),
        legend.title = element_blank()
  )

trust_pol <- ggplot(pol_outcomes %>% filter(outcome == "Trust"), 
                    aes(x=fct_rev(Vignette), y=mean, color = Pol, shape = Pol)) + 
  geom_point(size = 1.5) +
  geom_hline(yintercept=0, colour = "grey60", lty=2) + 
  scale_color_manual(values = c("#0015BC","#e9141d")) +
  scale_shape_manual(values = c(17,16) ) +
  coord_flip() +
  xlab("") + 
  ylab("\n") +
  labs(title="Trust") +
  ylim(-5,5) + 
  theme_bw() + 
  theme(plot.title = element_text(face="bold", hjust = 0.5),
        legend.position = "none",
        legend.background = element_rect(colour="grey80"),
        legend.title = element_blank()
  )

quality_pol <- ggplot(pol_outcomes %>% filter(outcome == "Quality"), 
                      aes(x=fct_rev(Vignette), y=mean, color = Pol, shape = Pol)) + 
  geom_point(size = 1.5) +
  geom_hline(yintercept=0, colour = "grey60", lty=2) + 
  scale_color_manual(values = c("#0015BC","#e9141d")) +
  scale_shape_manual(values = c(17,16) ) +
  coord_flip() +
  xlab("") + 
  ylab("\nMean Outcome") +
  labs(title="Quality") +
  ylim(-5,5) + 
  theme_bw() + 
  theme(plot.title = element_text(face="bold", hjust = 0.5),
        legend.position = "none",
        legend.background = element_rect(colour="grey80"),
        legend.title = element_blank()
  )

impact_pol <- ggplot(pol_outcomes %>% filter(outcome == "Impact"), 
                     aes(x=fct_rev(Vignette), y=mean, color = Pol, shape = Pol)) + 
  geom_point(size = 1.5) +
  geom_hline(yintercept=4, colour = "grey60", lty=2) + 
  scale_color_manual(values = c("#0015BC","#e9141d")) +
  scale_shape_manual(values = c(17,16) ) +
  coord_flip() +
  xlab("") + 
  ylab("\nMean Outcome") +
  labs(title="Impact") +
  ylim(2.5,5.5) + 
  theme_bw() + 
  theme(plot.title = element_text(face="bold", hjust = 0.5),
        legend.position = "none",
        legend.background = element_rect(colour="grey80"),
        legend.title = element_blank()
  )

jpeg("Figures/fig_5.jpg", units="in", width=9, height=7, res=300)
ggarrange(feeling_pol, trust_pol, quality_pol, impact_pol, nrow=2, ncol = 2, common.legend = T, legend = "bottom")
dev.off()
```

#Whether to Standardize Outcome Variables
```{r standardize variables}
#Main figures 3, 4, and 6 use standardized outcomes
standardize <- TRUE

#Because SDs and sample sizes are similar across treatment groups, we use Cohen's d
if(standardize == TRUE) {
  ads$Feeling <- scale(ads$Feeling)
  ads$Trust <- scale(ads$Trust)
  ads$Quality <- scale(ads$Quality)
  ads$Impact <- scale(ads$Impact)
}
```

####Reproduce Figure 3####
```{r figure 3}
#Perform regressions
out1_adjust <- lm(Feeling ~ Mechanism + Pol + Knowledge + Gender + Age + Race + Educ + Marital_Status + Income + Employment, data=ads)
out2_adjust <- lm(Trust ~ Mechanism + Pol + Knowledge + Gender + Age + Race + Educ + Marital_Status + Income + Employment, data=ads)
out3_adjust <- lm(Quality ~ Mechanism + Pol + Knowledge + Gender + Age + Race + Educ + Marital_Status + Income + Employment, data=ads)
out4_adjust <- lm(Impact ~ Mechanism + Pol + Knowledge + Gender + Age + Race + Educ + Marital_Status + Income + Employment, data=ads)

out1_adjust_tidy <- tidy(out1_adjust)[2:4,] %>% mutate(model = "Feeling")
out2_adjust_tidy <- tidy(out2_adjust)[2:4,] %>% mutate(model = "Trust")
out3_adjust_tidy <- tidy(out3_adjust)[2:4,] %>% mutate(model = "Quality")
out4_adjust_tidy <- tidy(out4_adjust)[2:4,] %>% mutate(model = "Impact")

main <- rbind(out1_adjust_tidy, out2_adjust_tidy, out3_adjust_tidy, out4_adjust_tidy) %>% arrange((model))
main$term <- rep(c("Bias","Lack of Transparency","Lack of Responsiveness"),4)
main$model <- factor(main$model, levels = c("Feeling","Trust","Quality","Impact"))
main <- main %>% arrange(model)

#Create figure
jpeg("Figures/fig_3.jpg", units="in", width=9, height=7, res=300)
dw_main <- dwplot(main %>% arrange(desc(model)),
vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2),
dot_args = list(aes(shape = model), size = 1.5)) +
scale_color_discrete(name = "Outcome", breaks=c("Feeling", "Trust", "Quality", "Impact")) +
  scale_shape_discrete(name = "Outcome", breaks=c("Feeling", "Trust", "Quality", "Impact")) +
    xlab("\nStandardized Treatment Effect") + ylab("") +
    xlim(-.75, .25) +
    theme_bw() + 
    theme(plot.title = element_text(face="bold", hjust = 0.5),
          legend.position = "right",
          legend.background = element_rect(colour="grey80")
          )

dw_main
dev.off()
```

#####Reproduce Figure 4####
```{r figure 4}
#Perform regressions
out1_adjust <- lm(Feeling ~ Vignette + Pol + Knowledge + Gender + Age + Race + Educ + Marital_Status + Income + Employment, data=ads)
out2_adjust <- lm(Trust ~ Vignette + Pol + Knowledge + Gender + Age + Race + Educ + Marital_Status + Income + Employment, data=ads)
out3_adjust <- lm(Quality ~ Vignette + Pol + Knowledge + Gender + Age + Race + Educ + Marital_Status + Income + Employment, data=ads)
out4_adjust <- lm(Impact ~ Vignette + Pol + Knowledge + Gender + Age + Race + Educ + Marital_Status + Income + Employment, data=ads)

out1_adjust_tidy <- tidy(out1_adjust)[2:4,] %>% mutate(model = "Child Welfare")
out2_adjust_tidy <- tidy(out2_adjust)[2:4,] %>% mutate(model = "Child Welfare")
out3_adjust_tidy <- tidy(out3_adjust)[2:4,] %>% mutate(model = "Child Welfare")
out4_adjust_tidy <- tidy(out4_adjust)[2:4,] %>% mutate(model = "Child Welfare")

out1_adjust_opp <- lm(Feeling ~ relevel(Vignette, ref="Courts Control") + Pol + Knowledge + Gender + Age + Race + Educ + Marital_Status + Income + Employment, data=ads)
out2_adjust_opp <- lm(Trust ~ relevel(Vignette, ref="Courts Control") + Pol + Knowledge + Gender + Age + Race + Educ + Marital_Status + Income + Employment, data=ads)
out3_adjust_opp <- lm(Quality ~ relevel(Vignette, ref="Courts Control") + Pol + Knowledge + Gender + Age + Race + Educ + Marital_Status + Income + Employment, data=ads)
out4_adjust_opp <- lm(Impact ~ relevel(Vignette, ref="Courts Control") + Pol + Knowledge + Gender + Age + Race + Educ + Marital_Status + Income + Employment, data=ads)
names(out1_adjust_opp$coefficients)[6:8] <- c("VignetteCourts Bias", "VignetteCourts Transparency", "VignetteCourts Responsiveness")
names(out2_adjust_opp$coefficients)[6:8] <- c("VignetteCourts Bias", "VignetteCourts Transparency", "VignetteCourts Responsiveness")
names(out3_adjust_opp$coefficients)[6:8] <- c("VignetteCourts Bias", "VignetteCourts Transparency", "VignetteCourts Responsiveness")
names(out4_adjust_opp$coefficients)[6:8] <- c("VignetteCourts Bias", "VignetteCourts Transparency", "VignetteCourts Responsiveness")

out1_adjust_opp_tidy <- tidy(out1_adjust_opp)[6:8,] %>% mutate(model = "Courts")
out2_adjust_opp_tidy <- tidy(out2_adjust_opp)[6:8,] %>% mutate(model = "Courts")
out3_adjust_opp_tidy <- tidy(out3_adjust_opp)[6:8,] %>% mutate(model = "Courts")
out4_adjust_opp_tidy <- tidy(out4_adjust_opp)[6:8,] %>% mutate(model = "Courts")

feeling <- rbind(out1_adjust_tidy, out1_adjust_opp_tidy) %>% arrange((model))
feeling$term <- rep(c("Bias","Lack of Transparency","Lack of Responsiveness"),2)

trust <- rbind(out2_adjust_tidy, out2_adjust_opp_tidy)
trust$term <- rep(c("Bias","Lack of Transparency","Lack of Responsiveness"),2)

quality <- rbind(out3_adjust_tidy, out3_adjust_opp_tidy)
quality$term <- rep(c("Bias","Lack of Transparency","Lack of Responsiveness"),2)

impact <- rbind(out4_adjust_tidy, out4_adjust_opp_tidy)
impact$term <- rep(c("Bias","Lack of Transparency","Lack of Responsiveness"),2)

#Create figure
jpeg("Figures/fig_4.jpg", units="in", width=9, height=7, res=300)

dw_feeling <- dwplot(arrange(feeling, desc(model)),
                     vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2),
                     dot_args = list(aes(shape = model)), size=1.5) +
  scale_color_manual(name = "model", breaks=c("Child Welfare","Courts"), values = c("#C77CFF","#7CAE00")) +
  scale_shape_discrete(name = "model", breaks=c("Child Welfare","Courts")) +
  xlab("\n ") + ylab("") +
  xlim(-.75, .25) +
  ggtitle("Feeling") +
  theme_bw() + 
  theme(plot.title = element_text(face="bold", hjust = 0.5),
        legend.justification=c(.6, 0),
        legend.background = element_rect(colour="grey80"),
        legend.title = element_blank()
  )

dw_trust <- dwplot(arrange(trust, desc(model)),
vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2),
                     dot_args = list(aes(shape = model)), size=1.5) +
  scale_color_manual(name = "model", breaks=c("Child Welfare","Courts"), values = c("#C77CFF","#7CAE00")) +
  scale_shape_discrete(name = "model", breaks=c("Child Welfare","Courts")) +
    xlab("\n ") + ylab("") +
    xlim(-.75, .25) +
    ggtitle("Trust") +
    theme_bw() + 
    theme(plot.title = element_text(face="bold", hjust = 0.5),
          legend.position = "none",
          legend.background = element_rect(colour="grey80"),
          legend.title = element_blank())

dw_quality <- dwplot(arrange(quality, desc(model)),
vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2),
                     dot_args = list(aes(shape = model)), size=1.5) +
  scale_color_manual(name = "model", breaks=c("Child Welfare","Courts"), values = c("#C77CFF","#7CAE00")) +
  scale_shape_discrete(name = "model", breaks=c("Child Welfare","Courts")) +
    xlab("\nStandardized Treatment Effect") + ylab("") +
    xlim(-.75, .25) +
    ggtitle("Quality") +
    theme_bw() + 
    theme(plot.title = element_text(face="bold", hjust = 0.5),
          legend.position = "none",
          legend.background = element_rect(colour="grey80"),
          legend.title = element_blank())

dw_impact <- dwplot(arrange(impact, desc(model)),
vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2),
                     dot_args = list(aes(shape = model)), size=1.5) +
  scale_color_manual(name = "model", breaks=c("Child Welfare","Courts"), values = c("#C77CFF","#7CAE00")) +
  scale_shape_discrete(name = "model", breaks=c("Child Welfare","Courts")) +
    xlab("\nStandardized Treatment Effect") + ylab("") +
    xlim(-.75, .25) +
    ggtitle("Impact") +
    theme_bw() + 
    theme(plot.title = element_text(face="bold", hjust = 0.5),
          legend.position = "none",
          legend.background = element_rect(colour="grey80"),
          legend.title = element_blank())

ggarrange(dw_feeling, dw_trust, dw_quality, dw_impact, nrow=2, ncol = 2, common.legend = T, legend = "bottom")

dev.off()
```

#####Z-tests to Evaluate Policy Sector Hypothesis####
```{r policy sector hypothesis} 
#Outcome 1
#Is the bias treatment effect for welfare significantly different from the bias treatment effect for courts?
z_score <- (as.numeric(coef(out1_adjust)[2]) - as.numeric(coef(out1_adjust_opp)[6]))/
  (sqrt((summary(out1_adjust)$coef[2,2])^2 + (summary(out1_adjust_opp)$coef[6,2])^2)) 
z.out1.bias <- 2*(1 - pnorm(abs(z_score)))

#Is the transparency treatment effect for welfare significantly different from the transparency treatment effect for courts?
z_score <- (as.numeric(coef(out1_adjust)[3]) - as.numeric(coef(out1_adjust_opp)[7]))/
  (sqrt((summary(out1_adjust)$coef[3,2])^2 + (summary(out1_adjust_opp)$coef[7,2])^2)) 
z.out1.trans <- 2*(1 - pnorm(abs(z_score)))

#Is the responsiveness treatment effect for welfare significantly different from the responsiveness treatment effect for courts?
z_score <- (as.numeric(coef(out1_adjust)[4]) - as.numeric(coef(out1_adjust_opp)[8]))/
  (sqrt((summary(out1_adjust)$coef[4,2])^2 + (summary(out1_adjust_opp)$coef[8,2])^2)) 
z.out1.responsiveness <- 2*(1 - pnorm(abs(z_score)))

#Outcome 2
#Is the bias treatment effect for welfare significantly different from the bias treatment effect for courts?
z_score <- (as.numeric(coef(out2_adjust)[2]) - as.numeric(coef(out2_adjust_opp)[6]))/
  (sqrt((summary(out2_adjust)$coef[2,2])^2 + (summary(out2_adjust_opp)$coef[6,2])^2)) 
z.out2.bias <- 2*(1 - pnorm(abs(z_score)))

#Is the transparency treatment effect for welfare significantly different from the transparency treatment effect for courts?
z_score <- (as.numeric(coef(out2_adjust)[3]) - as.numeric(coef(out2_adjust_opp)[7]))/
  (sqrt((summary(out2_adjust)$coef[3,2])^2 + (summary(out2_adjust_opp)$coef[7,2])^2)) 
z.out2.trans <- 2*(1 - pnorm(abs(z_score)))

#Is the responsiveness treatment effect for welfare significantly different from the responsiveness treatment effect for courts?
z_score <- (as.numeric(coef(out2_adjust)[4]) - as.numeric(coef(out2_adjust_opp)[8]))/
  (sqrt((summary(out2_adjust)$coef[4,2])^2 + (summary(out2_adjust_opp)$coef[8,2])^2)) 
z.out2.responsiveness <- 2*(1 - pnorm(abs(z_score)))

#Outcome 3
#Is the bias treatment effect for welfare significantly different from the bias treatment effect for courts?
z_score <- (as.numeric(coef(out3_adjust)[2]) - as.numeric(coef(out3_adjust_opp)[6]))/
  (sqrt((summary(out3_adjust)$coef[2,2])^2 + (summary(out3_adjust_opp)$coef[6,2])^2)) 
z.out3.bias <- 2*(1 - pnorm(abs(z_score)))

#Is the transparency treatment effect for welfare significantly different from the transparency treatment effect for courts?
z_score <- (as.numeric(coef(out3_adjust)[3]) - as.numeric(coef(out3_adjust_opp)[7]))/
  (sqrt((summary(out3_adjust)$coef[3,2])^2 + (summary(out3_adjust_opp)$coef[7,2])^2)) 
z.out3.trans <- 2*(1 - pnorm(abs(z_score)))

#Is the responsiveness treatment effect for welfare significantly different from the responsiveness treatment effect for courts?
z_score <- (as.numeric(coef(out3_adjust)[4]) - as.numeric(coef(out3_adjust_opp)[8]))/
  (sqrt((summary(out3_adjust)$coef[4,2])^2 + (summary(out3_adjust_opp)$coef[8,2])^2)) 
z.out3.responsiveness <- 2*(1 - pnorm(abs(z_score)))

#Outcome 4
#Is the bias treatment effect for welfare significantly different from the bias treatment effect for courts?
z_score <- (as.numeric(coef(out4_adjust)[2]) - as.numeric(coef(out4_adjust_opp)[6]))/
  (sqrt((summary(out4_adjust)$coef[2,2])^2 + (summary(out4_adjust_opp)$coef[6,2])^2)) 
z.out4.bias <- 2*(1 - pnorm(abs(z_score)))

#Is the transparency treatment effect for welfare significantly different from the transparency treatment effect for courts?
z_score <- (as.numeric(coef(out4_adjust)[3]) - as.numeric(coef(out4_adjust_opp)[7]))/
  (sqrt((summary(out4_adjust)$coef[3,2])^2 + (summary(out4_adjust_opp)$coef[7,2])^2)) 
z.out4.trans <- 2*(1 - pnorm(abs(z_score)))

#Is the responsiveness treatment effect for welfare significantly different from the responsiveness treatment effect for courts?
z_score <- (as.numeric(coef(out4_adjust)[4]) - as.numeric(coef(out4_adjust_opp)[8]))/
  (sqrt((summary(out4_adjust)$coef[4,2])^2 + (summary(out4_adjust_opp)$coef[8,2])^2)) 
z.out4.responsiveness <- 2*(1 - pnorm(abs(z_score)))

#Create table
z.policy.frame <- as.data.frame(matrix(c(z.out1.bias, z.out1.trans, z.out1.responsiveness,
                    z.out2.bias, z.out2.trans, z.out2.responsiveness,
                    z.out3.bias, z.out3.trans, z.out3.responsiveness,
                    z.out4.bias, z.out4.trans, z.out4.responsiveness), 3, 4))
row.names(z.policy.frame) <- c("Bias", "Transparency", "Responsiveness")
colnames(z.policy.frame) <- c("Feeling", "Trust", "Quality", "Impact")

#Results of z-tests
z.policy.frame 
```

####Reproduce Figure 6####
```{r conservative subgroup regressions}
#Perform regressions for conservatives
out1_adjust_conservative <- lm(Feeling ~ Vignette + Knowledge + Gender + Age + Race + Educ + Marital_Status + Income + Employment, data=ads %>% filter(Pol=="Conservative"))
out2_adjust_conservative <- lm(Trust ~ Vignette + Knowledge + Gender + Age + Race + Educ + Marital_Status + Income + Employment, data=ads %>% filter(Pol=="Conservative")) 
out3_adjust_conservative <- lm(Quality ~ Vignette + Knowledge + Gender + Age + Race + Educ + Marital_Status + Income + Employment, data=ads %>% filter(Pol=="Conservative"))
out4_adjust_conservative <- lm(Impact ~ Vignette + Knowledge + Gender + Age + Race + Educ + Marital_Status + Income + Employment, data=ads %>% filter(Pol=="Conservative"))

out1_adjust_conservative_tidy <- tidy(out1_adjust_conservative)[2:4,] %>% mutate(model = "Child Welfare")
out2_adjust_conservative_tidy <- tidy(out2_adjust_conservative)[2:4,] %>% mutate(model = "Child Welfare")
out3_adjust_conservative_tidy <- tidy(out3_adjust_conservative)[2:4,] %>% mutate(model = "Child Welfare")
out4_adjust_conservative_tidy <- tidy(out4_adjust_conservative)[2:4,] %>% mutate(model = "Child Welfare")

out1_adjust_conservative_opp <- lm(Feeling ~ relevel(Vignette, ref="Courts Control") + Knowledge + Gender + Age + Race + Educ + Marital_Status + Income + Employment, data=ads %>% filter(Pol=="Conservative"))
out2_adjust_conservative_opp <- lm(Trust ~ relevel(Vignette, ref="Courts Control") + Knowledge + Gender + Age + Race + Educ + Marital_Status + Income + Employment, data=ads %>% filter(Pol=="Conservative")) 
out3_adjust_conservative_opp <- lm(Quality ~ relevel(Vignette, ref="Courts Control") + Knowledge + Gender + Age + Race + Educ + Marital_Status + Income + Employment, data=ads %>% filter(Pol=="Conservative"))
out4_adjust_conservative_opp <- lm(Impact ~ relevel(Vignette, ref="Courts Control") + Knowledge + Gender + Age + Race + Educ + Marital_Status + Income + Employment, data=ads %>% filter(Pol=="Conservative"))

names(out1_adjust_conservative_opp$coefficients)[6:8] <- c("VignetteCourts Bias", "VignetteCourts Transparency", "VignetteCourts Responsiveness")
names(out2_adjust_conservative_opp$coefficients)[6:8] <- c("VignetteCourts Bias", "VignetteCourts Transparency", "VignetteCourts Responsiveness")
names(out3_adjust_conservative_opp$coefficients)[6:8] <- c("VignetteCourts Bias", "VignetteCourts Transparency", "VignetteCourts Responsiveness")
names(out4_adjust_conservative_opp$coefficients)[6:8] <- c("VignetteCourts Bias", "VignetteCourts Transparency", "VignetteCourts Responsiveness")

out1_adjust_conservative_opp_tidy <- tidy(out1_adjust_conservative_opp)[6:8,] %>% mutate(model = "Courts")
out2_adjust_conservative_opp_tidy <- tidy(out2_adjust_conservative_opp)[6:8,] %>% mutate(model = "Courts")
out3_adjust_conservative_opp_tidy <- tidy(out3_adjust_conservative_opp)[6:8,] %>% mutate(model = "Courts")
out4_adjust_conservative_opp_tidy <- tidy(out4_adjust_conservative_opp)[6:8,] %>% mutate(model = "Courts")

feeling_con <- rbind(out1_adjust_conservative_tidy, out1_adjust_conservative_opp_tidy) %>% arrange((model))
feeling_con$term <- rep(c("Bias","Lack of Transparency","Lack of Responsiveness"),2)

trust_con <- rbind(out2_adjust_conservative_tidy, out2_adjust_conservative_opp_tidy)
trust_con$term <- rep(c("Bias","Lack of Transparency","Lack of Responsiveness"),2)

quality_con <- rbind(out3_adjust_conservative_tidy, out3_adjust_conservative_opp_tidy)
quality_con$term <- rep(c("Bias","Lack of Transparency","Lack of Responsiveness"),2)

impact_con <- rbind(out4_adjust_conservative_tidy, out4_adjust_conservative_opp_tidy)
impact_con$term <- rep(c("Bias","Lack of Transparency","Lack of Responsiveness"),2)
```

```{r liberal subgroup regressions}
#Perform regressions for liberals
out1_adjust_liberal <- lm(Feeling ~ Vignette + Knowledge + Gender + Age + Race + Educ + Marital_Status + Income + Employment, data=ads %>% filter(Pol=="Liberal"))
out2_adjust_liberal <- lm(Trust ~ Vignette + Knowledge + Gender + Age + Race + Educ + Marital_Status + Income + Employment, data=ads %>% filter(Pol=="Liberal"))
out3_adjust_liberal <- lm(Quality ~ Vignette + Knowledge + Gender + Age + Race + Educ + Marital_Status + Income + Employment, data=ads %>% filter(Pol=="Liberal"))
out4_adjust_liberal <- lm(Impact ~ Vignette + Knowledge + Gender + Age + Race + Educ + Marital_Status + Income + Employment, data=ads %>% filter(Pol=="Liberal"))

out1_adjust_liberal_tidy <- tidy(out1_adjust_liberal)[2:4,] %>% mutate(model = "Child Welfare")
out2_adjust_liberal_tidy <- tidy(out2_adjust_liberal)[2:4,] %>% mutate(model = "Child Welfare")
out3_adjust_liberal_tidy <- tidy(out3_adjust_liberal)[2:4,] %>% mutate(model = "Child Welfare")
out4_adjust_liberal_tidy <- tidy(out4_adjust_liberal)[2:4,] %>% mutate(model = "Child Welfare")

out1_adjust_liberal_opp <- lm(Feeling ~ relevel(Vignette, ref="Courts Control") + Knowledge + Gender + Age + Race + Educ + Marital_Status + Income + Employment, data=ads %>% filter(Pol=="Liberal"))
out2_adjust_liberal_opp <- lm(Trust ~ relevel(Vignette, ref="Courts Control") + Knowledge + Gender + Age + Race + Educ + Marital_Status + Income + Employment, data=ads %>% filter(Pol=="Liberal"))
out3_adjust_liberal_opp <- lm(Quality ~ relevel(Vignette, ref="Courts Control") + Knowledge + Gender + Age + Race + Educ + Marital_Status + Income + Employment, data=ads %>% filter(Pol=="Liberal"))
out4_adjust_liberal_opp <- lm(Impact ~ relevel(Vignette, ref="Courts Control") + Knowledge + Gender + Age + Race + Educ + Marital_Status + Income + Employment, data=ads %>% filter(Pol=="Liberal"))

out1_adjust_liberal_opp_tidy <- tidy(out1_adjust_liberal_opp)[6:8,] %>% mutate(model = "Courts")
out2_adjust_liberal_opp_tidy <- tidy(out2_adjust_liberal_opp)[6:8,] %>% mutate(model = "Courts")
out3_adjust_liberal_opp_tidy <- tidy(out3_adjust_liberal_opp)[6:8,] %>% mutate(model = "Courts")
out4_adjust_liberal_opp_tidy <- tidy(out4_adjust_liberal_opp)[6:8,] %>% mutate(model = "Courts")

feeling_lib <- rbind(out1_adjust_liberal_tidy, out1_adjust_liberal_opp_tidy) %>% arrange((model))
feeling_lib$term <- rep(c("Bias","Lack of Transparency","Lack of Responsiveness"),2)

trust_lib <- rbind(out2_adjust_liberal_tidy, out2_adjust_liberal_opp_tidy)
trust_lib$term <- rep(c("Bias","Lack of Transparency","Lack of Responsiveness"),2)

quality_lib <- rbind(out3_adjust_liberal_tidy, out3_adjust_liberal_opp_tidy)
quality_lib$term <- rep(c("Bias","Lack of Transparency","Lack of Responsiveness"),2)

impact_lib <- rbind(out4_adjust_liberal_tidy, out4_adjust_liberal_opp_tidy)
impact_lib$term <- rep(c("Bias","Lack of Transparency","Lack of Responsiveness"),2)

names(out1_adjust_liberal_opp$coefficients)[6:8] <- c("VignetteCourts Bias", "VignetteCourts Transparency", "VignetteCourts Responsiveness")
names(out2_adjust_liberal_opp$coefficients)[6:8] <- c("VignetteCourts Bias", "VignetteCourts Transparency", "VignetteCourts Responsiveness")
names(out3_adjust_liberal_opp$coefficients)[6:8] <- c("VignetteCourts Bias", "VignetteCourts Transparency", "VignetteCourts Responsiveness")
names(out4_adjust_liberal_opp$coefficients)[6:8] <- c("VignetteCourts Bias", "VignetteCourts Transparency", "VignetteCourts Responsiveness")
```

```{r figure 6}
#Create figure
feeling <- rbind(feeling_con, feeling_lib) %>% mutate(pol = rep(c("Conservative","Liberal"), each = 6)) %>% 
  mutate(model = paste0(model," ",pol))
feeling$model <- factor(feeling$model, levels = c("Child Welfare Conservative","Courts Conservative","Child Welfare Liberal","Courts Liberal"))

trust <- rbind(trust_con, trust_lib) %>% mutate(pol = rep(c("Conservative","Liberal"), each = 6)) %>% 
  mutate(model = paste0(model," ",pol))
trust$model <- factor(trust$model, levels = c("Child Welfare Conservative","Courts Conservative","Child Welfare Liberal","Courts Liberal"))

quality <- rbind(quality_con, quality_lib) %>% mutate(pol = rep(c("Conservative","Liberal"), each = 6)) %>% 
  mutate(model = paste0(model," ",pol))
quality$model <- factor(quality$model, levels = c("Child Welfare Conservative","Courts Conservative","Child Welfare Liberal","Courts Liberal"))

impact <- rbind(impact_con, impact_lib) %>% mutate(pol = rep(c("Conservative","Liberal"), each = 6)) %>% 
  mutate(model = paste0(model," ",pol))
impact$model <- factor(impact$model, levels = c("Child Welfare Conservative","Courts Conservative","Child Welfare Liberal","Courts Liberal"))

jpeg("Figures/fig_6.jpg", units="in", width=9, height=7, res=300)

combined_feeling <- dwplot(arrange(feeling, desc(model)),
                      vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2),
                      dot_args = list(aes(shape = model, color = model)), size=1.5) +
  scale_color_manual(name = "model", breaks=c("Child Welfare Conservative","Courts Conservative","Child Welfare Liberal","Courts Liberal"), values = c("firebrick4","firebrick3","dodgerblue4","dodgerblue1")) +
  scale_shape_manual(name = "model", breaks=c("Child Welfare Conservative","Courts Conservative","Child Welfare Liberal","Courts Liberal"), values = c(17,16,17,16)) +
  xlab("\n ") + ylab("") +
  xlim(-1.0, .5) +
  ggtitle("Feeling") +
  theme_bw() + 
  theme(plot.title = element_text(face="bold", hjust = 0.5),
        legend.justification=c(.6, 0),
        legend.background = element_rect(colour="grey80"),
        legend.title = element_blank()
  )

combined_trust <- dwplot(arrange(trust, desc(model)),
                      vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2),
                      dot_args = list(aes(shape = model, color = model)), size=1.5) +
  scale_color_manual(name = "model", breaks=c("Child Welfare Conservative","Courts Conservative","Child Welfare Liberal","Courts Liberal"), values = c("firebrick4","firebrick3","dodgerblue4","dodgerblue1")) +
  scale_shape_manual(name = "model", breaks=c("Child Welfare Conservative","Courts Conservative","Child Welfare Liberal","Courts Liberal"), values = c(17,16,17,16)) +
  xlab("\n ") + ylab("") +
  xlim(-1.0, .5) +
  ggtitle("Trust") +
  theme_bw() + 
  theme(plot.title = element_text(face="bold", hjust = 0.5),
        #legend.justification=c(.6, 0),
        legend.position = "none",
        legend.background = element_rect(colour="grey80"),
        legend.title = element_blank()
  )

combined_quality <- dwplot(arrange(quality, desc(model)),
                      vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2),
                      dot_args = list(aes(shape = model, color = model)), size=1.5) +
  scale_color_manual(name = "model", breaks=c("Child Welfare Conservative","Courts Conservative","Child Welfare Liberal","Courts Liberal"), values = c("firebrick4","firebrick3","dodgerblue4","dodgerblue1")) +
  scale_shape_manual(name = "model", breaks=c("Child Welfare Conservative","Courts Conservative","Child Welfare Liberal","Courts Liberal"), values = c(17,16,17,16)) +
  xlab("\nStandardized Treatment Effect") + ylab("") +
  xlim(-1.0, .5) +
  ggtitle("Quality") +
  theme_bw() + 
  theme(plot.title = element_text(face="bold", hjust = 0.5),
        #legend.justification=c(.6, 0),
        legend.position = "none",
        legend.background = element_rect(colour="grey80"),
        legend.title = element_blank()
  )

combined_impact <- dwplot(arrange(impact, desc(model)),
                      vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2),
                      dot_args = list(aes(shape = model, color = model)), size=1.5) +
  scale_color_manual(name = "model", breaks=c("Child Welfare Conservative","Courts Conservative","Child Welfare Liberal","Courts Liberal"), values = c("firebrick4","firebrick3","dodgerblue4","dodgerblue1")) +
  scale_shape_manual(name = "model", breaks=c("Child Welfare Conservative","Courts Conservative","Child Welfare Liberal","Courts Liberal"), values = c(17,16,17,16)) +
  xlab("\nStandardized Treatment Effect") + ylab("") +
  xlim(-1.0, .5) +
  ggtitle("Impact") +
  theme_bw() + 
  theme(plot.title = element_text(face="bold", hjust = 0.5),
        #legend.justification=c(.6, 0),
        legend.position = "none",
        legend.background = element_rect(colour="grey80"),
        legend.title = element_blank()
  )

ggarrange(combined_feeling, combined_trust, combined_quality, combined_impact, nrow=2, ncol = 2, common.legend = T, legend = "bottom")

dev.off()
```

#####Z-tests to Evaluate Political Ideology Hypothesis####
```{r political ideology hypothesis}
#Outcome 1
#Is the courts bias treatment effect for liberals significantly different from the courts bias treatment effect for conservatives?
z_score <- (as.numeric(coef(out1_adjust_liberal_opp)[2]) - as.numeric(coef(out1_adjust_conservative_opp)[6]))/
  (sqrt((summary(out1_adjust_liberal_opp)$coef[2,2])^2 + (summary(out1_adjust_conservative_opp)$coef[6,2])^2)) 
z.out1.bias <- 2*(1 - pnorm(abs(z_score)))

#Is the courts transparency treatment effect for liberals significantly different from the courts transparency treatment effect for conservatives?
z_score <- (as.numeric(coef(out1_adjust_liberal_opp)[3]) - as.numeric(coef(out1_adjust_conservative_opp)[7]))/
  (sqrt((summary(out1_adjust_liberal_opp)$coef[3,2])^2 + (summary(out1_adjust_conservative_opp)$coef[7,2])^2)) 
z.out1.trans <- 2*(1 - pnorm(abs(z_score)))

#Is the courts responsiveness treatment effect for liberals significantly different from the courts responsiveness treatment effect for conservatives?
z_score <- (as.numeric(coef(out1_adjust_liberal_opp)[4]) - as.numeric(coef(out1_adjust_conservative_opp)[8]))/
  (sqrt((summary(out1_adjust_liberal_opp)$coef[4,2])^2 + (summary(out1_adjust_conservative_opp)$coef[8,2])^2)) 
z.out1.responsiveness <- 2*(1 - pnorm(abs(z_score)))

#Outcome 2
#Is the courts bias treatment effect for liberals significantly different from the courts bias treatment effect for conservatives?
z_score <- (as.numeric(coef(out2_adjust_liberal_opp)[2]) - as.numeric(coef(out2_adjust_conservative_opp)[6]))/
  (sqrt((summary(out2_adjust_liberal_opp)$coef[2,2])^2 + (summary(out2_adjust_conservative_opp)$coef[6,2])^2)) 
z.out2.bias <- 2*(1 - pnorm(abs(z_score)))

#Is the courts transparency treatment effect for liberals significantly different from the courts transparency treatment effect for conservatives?
z_score <- (as.numeric(coef(out2_adjust_liberal_opp)[3]) - as.numeric(coef(out2_adjust_conservative_opp)[7]))/
  (sqrt((summary(out2_adjust_liberal_opp)$coef[3,2])^2 + (summary(out2_adjust_conservative_opp)$coef[7,2])^2)) 
z.out2.trans <- 2*(1 - pnorm(abs(z_score)))

#Is the courts responsiveness treatment effect for liberals significantly different from the courts responsiveness treatment effect for conservatives?
z_score <- (as.numeric(coef(out2_adjust_liberal_opp)[4]) - as.numeric(coef(out2_adjust_conservative_opp)[8]))/
  (sqrt((summary(out2_adjust_liberal_opp)$coef[4,2])^2 + (summary(out2_adjust_conservative_opp)$coef[8,2])^2)) 
z.out2.responsiveness <- 2*(1 - pnorm(abs(z_score)))

#Outcome 3
#Is the courts bias treatment effect for liberals significantly different from the courts bias treatment effect for conservatives?
z_score <- (as.numeric(coef(out3_adjust_liberal_opp)[2]) - as.numeric(coef(out3_adjust_conservative_opp)[6]))/
  (sqrt((summary(out3_adjust_liberal_opp)$coef[2,2])^2 + (summary(out3_adjust_conservative_opp)$coef[6,2])^2)) 
z.out3.bias <- 2*(1 - pnorm(abs(z_score)))

#Is the courts transparency treatment effect for liberals significantly different from the courts transparency treatment effect for conservatives?
z_score <- (as.numeric(coef(out3_adjust_liberal_opp)[3]) - as.numeric(coef(out3_adjust_conservative_opp)[7]))/
  (sqrt((summary(out3_adjust_liberal_opp)$coef[3,2])^2 + (summary(out3_adjust_conservative_opp)$coef[7,2])^2)) 
z.out3.trans <- 2*(1 - pnorm(abs(z_score)))

#Is the courts responsiveness treatment effect for liberals significantly different from the courts responsiveness treatment effect for conservatives?
z_score <- (as.numeric(coef(out3_adjust_liberal_opp)[4]) - as.numeric(coef(out3_adjust_conservative_opp)[8]))/
  (sqrt((summary(out3_adjust_liberal_opp)$coef[4,2])^2 + (summary(out3_adjust_conservative_opp)$coef[8,2])^2)) 
z.out3.responsiveness <- 2*(1 - pnorm(abs(z_score)))

#Outcome 4
#Is the courts bias treatment effect for liberals significantly different from the courts bias treatment effect for conservatives?
z_score <- (as.numeric(coef(out4_adjust_liberal_opp)[2]) - as.numeric(coef(out4_adjust_conservative_opp)[6]))/
  (sqrt((summary(out4_adjust_liberal_opp)$coef[2,2])^2 + (summary(out4_adjust_conservative_opp)$coef[6,2])^2)) 
z.out4.bias <- 2*(1 - pnorm(abs(z_score)))

#Is the courts transparency treatment effect for liberals significantly different from the courts transparency treatment effect for conservatives?
z_score <- (as.numeric(coef(out4_adjust_liberal_opp)[3]) - as.numeric(coef(out4_adjust_conservative_opp)[7]))/
  (sqrt((summary(out4_adjust_liberal_opp)$coef[3,2])^2 + (summary(out4_adjust_conservative_opp)$coef[7,2])^2)) 
z.out4.trans <- 2*(1 - pnorm(abs(z_score)))

#Is the courts responsiveness treatment effect for liberals significantly different from the courts responsiveness treatment effect for conservatives?
z_score <- (as.numeric(coef(out4_adjust_liberal_opp)[4]) - as.numeric(coef(out4_adjust_conservative_opp)[8]))/
  (sqrt((summary(out4_adjust_liberal_opp)$coef[4,2])^2 + (summary(out4_adjust_conservative_opp)$coef[8,2])^2)) 
z.out4.responsiveness <- 2*(1 - pnorm(abs(z_score)))

#Create table
z.subgroup.frame <- as.data.frame(matrix(c(z.out1.bias, z.out1.trans, z.out1.responsiveness,
                    z.out2.bias, z.out2.trans, z.out2.responsiveness,
                    z.out3.bias, z.out3.trans, z.out3.responsiveness,
                    z.out4.bias, z.out4.trans, z.out4.responsiveness), 3, 4))
row.names(z.subgroup.frame) <- c("Bias", "Transparency", "Responsiveness")
colnames(z.subgroup.frame) <- c("Feeling", "Trust", "Quality", "Impact")

#Results of z-tests
z.subgroup.frame
```

